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Conventional methods for the simulation of diffusive systems are quite slow when applied to 
strongly inhomogeneous systems. We present a new hierarchical approach based on dynamic 
renormalization-group ideas and on the Walsh transform (or Haar wavelet) of signal-processing 
theory. The method is very efficient for simulation of petroleum reservoirs or other strongly inho- 
mogeneous diffusive or pressure-driven flow systems. In a test case, the hierarchical method is found 
to achieve 1.5% accuracy roughly 25 times faster than conventional finite difference methods. 

so- 
on ; 

0^ | Traditional finitc-clcmcnt and finite-difference methods for numerical solution of differential equations have a dis- 
cretization error that depends on a power of the time or space increment, At or Ar. In the case of a diffusive system, 
O ' stability usually requires At to be of order Ar 2 /D (where D is the diffusivity) , and the discretization error is propor- 
tional to a power of Ar. In highly inhomogeneous systems, parts of which require very small Ar and/or very large 
' diffusivity D, this requires a very small At. 

The instability for large At arises when material moves more than one cell diameter during the time interval At; 
£N) ' in the usual explicit finite difference (EFD) algorithm, material is allowed to move from a cell only into its nearest 
neighbors. Our approach solves this problem by allowing material to move across as many cells as necessary. We 

■ describe the motion of the material by a discrete Green function or influence function G\t{d, s) such that GAt(d, s)c(s) 
is the amount of material that moves from a source cell s (whose original material content is c(s)) to a destination cell 
d during the time interval At. If GAt(d, s) is nonzero only when d and s are nearest neighbor cells, this is equivalent 
to a conventional EFD algorithm. This is the case for sufficiently small At, so we may begin with such an algorithm 

, and coarsen the time scale by doubling At. The Green function for the interval 2At is obtained from that for At by 

■ self-convolution (see Eq. |^ below) . 

SO \ Of course, this repeated convolution process increases the spatial range of the influence function, which is stored in 
our computer implementation as a linked list; soon the number of destination cells d that can be reached from each 
source cell s becomes large and the method becomes very time-consuming. However, we can coarsen the space as well 
as the time scale, by lumping cells together into larger cells. This decreases the number of source cells, as well as the 
number of destination cells reached from each source cell, and hence the size of the Green- function list. This scheme 
is motivated by the renormalization-group method which has been so useful in the theory of critical phenomena [[l], 
although our present description does not require prior knowledge of renormalization-group theory. It has been shown 
[p| that the diffusion problem in a homogeneous system has a fixed point with respect to a combined space-and-time 
renormalization transformation. That is, we can continue coarsening the space and time scales indefinitely without 
indefinitely increasing the size of the Green-function list. 

Such cell coarsening is even more useful in an inhomogeneous system, because we can use physical information to 
choose which cells to lump together. An inhomogeneous system such as an oil reservoir tends to be compartmentalized 
into compartments within which oil flows fairly freely, separated by relatively impermeable regions. If we lump cells 
between which there is relatively free flow (high effective diffusivity), when we reach a large scale the cells will be the 
compartments. 

We can visualize the hierarchical lumping of cells by placing the cells on a binary tree as in Fig. El 
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FIG. 1. Sketch of a hierarchically subdivided system (left) and its representation as a binary tree (right). Cell L is subdivided 
into I and I' as described in the text. 



The disadvantage of this coarsening process, of course, is that we lose spatial resolution in our description of the 
system. To some extent such a contraction of the description is desired, but we would like to maintain control of 
our approximations and limit the loss of information. When we replace contents c{l) and c(/'), of two cells I and V 
by the coarse-grained content c(L) = c(l) + c(V) where the larger cell L is the union of I and we can avoid losing 
any information if we also include as a variable the difference c(L,l) = c{l) — c(V) as well as the sum c(L,0) (the 
arguments and 1 simply indicate whether a sum or a difference is intended). 

We can do this at any level of the tree shown in Fig. |l| - the difference between the two halves of the small cell 
I can be denoted by c(l, 1). We can even define a difference of differences c(L, 1, 1) = c(l, 1) — c{V , 1). Each "1" in 
this expression can be regarded as one bit of a "Walsh sequency index" w. (The Walsh transform is a discrete signal 
transform used in electrical engineering || which we here generalize to a hierarchical system.) We will lump the 
bits into a single binary integer (Walsh index) u>, so c(L, 1,1) becomes c(L,w) with w — 3 (11 in binary). We can 
define a general "Walsh variable" recursively by 

c(L, wO) = c(/, w) + c(l', w) 

c(L,wl)=c(l,w)-c(l',w) ^> 

whenever I and I' are the two halves of the cell L. Here the Walsh index wO is w with a zero bit appended at 
the right, i.e., wO = 2w; similarly wl = 2w + 1. In this notation, our original c(l) is written c(l,0) and serves to 
start the recursion. The Walsh index plays a role similar to that of the wavenumber in the Fourier transform, in 
that variables with a small Walsh index describe large-scale structure, whereas large Walsh indices describe small 
wavelength structure within a cell L. For each "1" bit in the binary representation of w, there is one subtraction in 
the construction of c(l, w). 

If each cell-lumping replaces two cell contents by two Walsh variables, we never lose variables and our equations 
remain as complicated as ever (but exact). However, we are now in a position to selectively drop small terms in the 
Green function: the variables with large Walsh indices (i.e., many differences rather than sums) will be less important 
than those with small Walsh indices. Our algorithm drops terms from the Green function if they are less than some 
preset tolerance 8. Typically, terms involving differences are much smaller than those involving sums, so these are 
the ones dropped. This is the virtue of the hierarchical description: terms describing the effects of a cell content are 
never negligible compared to the terms for nearby cells, whereas terms describing the effects of differences may be 
negligible compared to terms for sums. These advantages are similar to those of spectral (Fourier transform) methods 
in homogeneous systems; in a sense one can regard the Walsh transform as the proper generalization of the Fourier 
transform to inhomogeneous systems. 

Commercial reservoir simulation programs usually treat a reservoir as a three-component system (oil, gas, 
water) in which flow is governed by Darcy's law. However, to provide a straightforward test of the hierarchical 
method described above, we will consider only one component (oil). In this case, Darcy's law reduces to the diffusion 
equation 

^M=V-[D(r)VP(r,t)] (2) 

where the pressure P plays the role of the diffusing density; the effective diffusivity D{y) is proportional to the perme- 
ability. So the problem we will actually solve is that of diffusion in a very inhomogeneous system; the inhomogeneity 
is contained in the function D(r). 
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To describe the evolution of the cell content Ct(d) of a cell labeled d (proportional to the "density" P), the most 
straightforward discretization of Eq. @ is 



c t +At{d) - c t (d) 



= Ar-2j2 D (f)lct(d + (f))-c t (d)} (3) 



where the sum is over faces / of the cell d, and d+ (/) is the cell in front of the (directed) face / (the cell behind it is 
always d). When we lump cells, so some of our variables are Walsh variables c(d,w), we can write the discretization 
(Eq. ||) in the form 

Ct+At(d,w) = ^2 G At (d,w; s,v)c t (s,v) (4) 

where GAt(d, w; s, v) is a Green function or influence function describing the influence of a Walsh variable in the 
source cell s on one in the destination cell d. Before any lumping has occurred, all Walsh variables have w = 0, and 
GAt(d, w; s, v) = (At/ Ar 2 )D(f) if d and s are neighbors separated by the face / and w — v — 0, and zero otherwise. 
Equation His thus equivalent to an EFD algorithm, which requires a small At; we have used the practical limit of 
stability M, 

At = Ar 2 /4A„ ax (5) 

where D max is the maximum diffusivity in the system. 

We now increase At to 2At; the new Green function is the convolution 

G 2 At(d, w; 8,v) = y~] GAt{d, w; e, u)GAt(e, u; s, v) (6) 

e,u 

After several such convolutions, the spatial range of the Green function is increased, especially in regions of high 
diffusivity. Here the pressure (i.e., density) in nearby cells equalizes rapidly - the contents c(l) and c(V) of two nearby 
cells contribute nearly equally to future contents c(d): G(d, 0;/,0) « G(d, 0;Z',0) (the zeroes here are the Walsh 
indices). To decide whether two cells / and I' should be lumped together, we look at the ratio (using d = I) 

G(Z,0;Z',0) 



Fig. Jj) 



we can calculate the new Green function 



G(Z,0;Z,0) 

Generally r < 1; we have used r > 0.91 as our lumping criterion (see Fig 

When we decide to lump two cells I and I' into a large cell L (as in Fig 
in two stages. In the first stage we calculate elements G'(ci, w; s, v) in which the destination cell d takes coarse values 
(including L) but s takes values including I and V . These are the same as the old Gs unless d is L, in which case we 
obtain from Eq. |l| 

G'(L, wO; s, v) = G(l, w; s, v) + G(V , w; s, v) 

G'(L, wl; s, v) = G(l, w; s, v) - G(l', w; s, v) W 

where as before wO means w with a zero bit appended. In the second stage, we calculate G"(d, w; s, v) where both d 
and s take values L and not I or V : we coarsen the source cell. Again, G" (d, w; s, v) — G'(d, w] s, v) unless s — L, in 
which case 



G"(d,w;L,vQ) = ±[G'(d,w;s,v) + G'(d,w;s',v)} 
G"(d,w;L,vl) = ±[G'(d,w;s,v) - G'(d,w;s' ,v)] 



(9) 



Although we have implemented our algorithm in three dimensions, we use a 2D system for our test calculation The 
test system has four control parameters: N, I, 8, and r. The first two describe the complexity of the system: N is 
the system size (16 x 16, 32 x 32, or 64 x 64) and the inhomogeneity parameter J is a normalized standard deviation: 
the standard deviation of the permeability divided by the mean permeability. The other two parameters are the error 
tolerance S, which we choose to give an acceptable overall truncation error, and the lumping threshold r (Eq. [?]), 
which we tune to maximize the speed of the algorithm. 

The diffusivity (i.e., permeability) distribution we have used is a realization of a log-normal distribution with fractal 
spatial correlations, obtained by exponentiating a correlated Gaussian distribution described elsewhere 0. We adjust 
the normalized standard deviation I by scaling the Gaussian distribution before exponentiating it. The prefactor that 



3 



governs the overall scale of the diffusivity or permeability can be removed from the problem by rescaling time. The 
system shown in Fig. ^|is 64 x 64; we specify the permeability in the smaller- N systems by coarse- graining (averaging 
over 2 x 2 or 4 x 4 cells). The permeability at a face is taken to be the average of that in the adjoining cells. 

As a test initial condition, we use a delta- function density concentrated in the lower left cell of the system. After 
an infinite time, the density takes a uniform value P^; we evolve the system until the density in the source cell is 
2Poo, as shown in Fig. ||. 
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FIG. 2. The final density distribution in the 64 x 64 system with inhomogeneity I — 20, showing boundaries of lumped cells. 

To compare our scheme with an EFD algorithm, we look first at the 64 x 64 system with a substantial inhomogeneity 
(/ = 20), and use r = 0.9 (justified by Fig. ^ below). We vary the remaining parameter, the tolerance 5, and plot 
in Fig. I) the required CPU time (on a Silicon Graphics R4000PC Indy, 133 MHz) against the accuracy achieved, 
defined by the fractional rms error 

error 2 = A^ 1 ^[ P(c) ~ Pexact{c) f. (10) 
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FIG. 3. Logarithmic plot of the CPU time required by our hierarchical Green function algorithm, compared to that required 
by an explicit finite difference (EFD) algorithm, for a 64 x 64 system with inhomogeneity / = 20. The rightmost data point for 
each system size has tolerance parameter S — 10~ 2 , and those for 32 x 32 system decrease by factors of 10. For 64 x 64 there 
are points at 5 x 10~ n as well as at 10~". 

Note that the speedup factor of our algorithm compared to the EFD algorithm increases rapidly as the allowable 
error is increased. It is indicated by a double arrow at the error value of 1.5%, where it is about 25. Using this 
factor as a figure of merit for our algorithm, we plot it in Fig. |as a function of the lumping parameter r; evidently 
performance improves as we increase r toward 1.0. Placing r very close to 1.0 risks magnifying the effects of small 
numerical errors, so we used r = 0.9 in the other figures. 
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FIG. 4. Speedup factor as a function of the lumping threshold r (Eq. Q), in a 32 x 32 system with inhomogeneity / = 20 
and tolerance S = 10~ 4 . 
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We are using the explicit finite difference algorithm as our point of comparison, not because it is the most efficient 
existing algorithm, but because it is the simplest. However, the use of more efficient implicit algorithms increases the 
allowed At only by factors of order 2 and does not affect our conclusions. 

The dependence of the speedup factor on the system properties is shown in Table |. Evidently the Green function 
algorithm is most advantageous in highly inhomogeneous systems. 
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TABLE I. Speedup factor as a function of system size N, for three values of the inhomogeneity /. Tolerance is S = 10 



Although the method described here has some features in common with methods already in common use in grid- 
based numerical simulation, none of these older methods approaches its efficiency for inhomogeneous systems. The 
idea of coarsening cells is used in "multigrid" methods ||. For example, the solution of Laplace's equation by the 
relaxation method is very slow on a fine grid, ft can be speeded up by doing a few iterations of relaxation on a larger 
grid to get the coarse features of the solution correct, and then returning to the fine grid to improve the finer features. 
Computational fluid dynamics codes often use an "adaptive grid" method ||, wherein larger grid sizes are used in 
regions where fields do not vary rapidly in space, and finer grids are used where there are fine-scale variations in the 
fields. These regions are typically rectangular and cannot follow compartment shapes as closely as ours. Unlike in 
our approach, the time increment cannot be increased above what is stable on the finest grid. 

In conclusion, we have shown that a hierarchical algorithm based on the dynamic renormalization group and the 
Walsh transform can simulate diffusive flow in an inhomogeneous system much more efficiently than conventional 
finite-difference algorithms. This occurs because the rapid power-law dependence of CPU time on the fundamental 
scales Ar and At is replaced in the hierarchical method by a logarithmic dependence. 
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